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Abs tract 


Summer research at NASA Lewis Research Center gave the opportunity 
to incorporate new control volumes in the Denton 3-D finite-volume time- 
marching code. For duct flows, the new control volumes require no 
transverse smoothing and this allows calculations with large transverse 
gradients in properties without significant numerical total pressure 
losses. 

The summer research also pointed to possibilities for improving the 
Denton code to obtain better distributions of properties through 
shocks. Much better total pressure distributions through shocks are 
obtained when the interpolated effective pressure, needed to stabilize 
the solution procedure, is used to calculate the total pressure. This 
simple change largely eliminates the undershoot in total pressure down- 
stream of a shock. Overshoots and undershoots in total pressure can 
then be further reduced by a factor of 10 by adopting the effective 
density method, developed at VPI&SU, rather than the effective pressure 
method. Use of a Mach number dependent interpolation scheme for pres- 
sure then removes the overshoot in static pressure downstream of a 
shock. 

The stability of interpolation schemes used for the calculation of 
effective density is analyzed and a Mach number dependent scheme, the 
M&M formula, is developed. This formula combines the advantages of the 
correct perfect gas equation for subsonic flow with the stability of 2- 
polnt and 3-point interpolation schemes for supersonic flow. 
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PART 1 


SUMMER WORK AT NASA LEWIS RESEARCH CENTER 


INVISCID FLOW CALCULATIONS USING THE DENTON CODE 


1.1 NEW CONTROL VOLUMES IN THREE DIMENSIONS 

A new control volume has been introduced (1) which allows the 
calculation of transonic flow in ducts using the f inite- volume 
method without the smoothing of flow properties that is usually 
needed(2). Previous work(l) using these new control volumes has 
been limited to two dimensions. The first part of the work at 
NASA Lewis this summer was to extend these new control volumes to 
three-dimensional flow calculations. This was thought important 
since the three-dimensional version of the Denton f inite- volume 
code is the one typically used at NASA. 

An example of a typical new three dimensional control volume 
is shown in Fig. 1. The locations of control volume boundaries 
are specified in the input data and the control volume surfaces 
are constructed from this information. Once the control volume 
boundaries are known then the grid points are placed in the 
middle of the upstream and downstream faces of the control 
volume. The fluxes through the transverse faces of the control 
volume needed for the continuity and momentum balances are 
determined from interpolated properties using the nodes adjacent 
to the face. Fig. 2 shows two adjacent control volumes of 
different sizes. The procedure for calculating the properties to 
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be used in calculating the fluxes on the common boundary ( face I) 
can be shown in the following way. For face I, any property, X, 
is determined from the average of the property at points A and 
B, where the values of the properties X A and X fi are determined by 
linearly interpolating between the values of the property at 
nodes 1 and 2 , and between the values of the property at nodes 3 
and 4, respectively. 

Assuming that face II corresponds to a solid boundary, the 
values of a property at points C and D are determined by linear 
extrapolation using the values of the property at nodes 1 and 2 , 
and 3 and 4 ; respectively. For the present calculations, only 
the pressure needs to be calculated at the solid boundaries since 
the fluxes of mass are set equal to zero through these solid 
boundaries. 

Additional adjustments were made to NASA's finite volume 
code to allow the calculation of cascade geometries with the new 
control volumes . Fig. 3 shows a two dimensional projection of a 
typical grid system up to the leading edge of a cascade blade. 
Note that a grid point is not located along the periodic boundary 
when the new control volumes are used . The computational domain 
extends from the lower periodic boundary to the upper periodic 
boundary. The missing calculation points outside the 
computational domain are replaced by the corresponding points 
adjacent to the other periodic boundary. For the calculations 
made this summer the leading and trailing edges were modeled as 
shown in Fig. 4. 
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1.2 SAJBEN'S DIFFUSER CALCULATIONS USING THE NEW CONTROL VO IDMES 


Sajben's diffuser(3) was used as a test case to test the 
effect which the new control volumes have on the calculation of 
transonic flow. The results from calculations using NASA's 
current f inite- volume code are used for comparison. 

The geometry and grid used in the calculations are shown in 
Fig. 5. There were 34 axial grid points and 10 equally spaced 
radial grid points. The current NASA code requires input in 
x-r-0 coordinates. This requires that the two dimensional nozzle 
geometry be input either in the x-0 coordinates or in the x-r 
coordinates. Both were used successfully. The current 
calculations are made essentially two dimensional by inputing the 

coordinates of the diffuser at a very large radius (900 m.) in 

{ 

x-r coordinates. The calculations begin at x/h=-3.6 and end at 
x/h=7.9, where h is the throat height. The inlet total pressure 
is 135 kPa and the inlet total temperature is 300 K. The exit 
static pressure is 108 kPa. The gives a p ex i t / pt i n i e t~° • 800 • 

With these conditions, one dimensional isentropic flow gives a 
shock with an upstream Mach number of 1.495 at the location 
marked in Fig. 5. Multigridding is used to improve the 
convergence speed. A copy of the input file used for these 
calculations is in Appendix A. 

■ /• 

Fig. 6 shows a comparison of bottom flat wall static 
pressures obtained using the old control volumes and using the 
new control volumes( "old" will refer to the type of control 
volumes used in NASA’s current finite volume code and "new" will 
refer to the type of control volumes shown in Fig. 1). A one 
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dimensional analytical calculation was made for the Sajben 
diffuser geometry using the above specified boundary conditions. 
Fig. 6 includes the static pressure distribution of this 
one- dimensional calculation for comparison with the numerical 
solutions. The agreement is good except through the shock. 

1.3 EV AID ATI ON OF TOTAL PRESSURE IN THE DENTON CODE 

Fig. 7 shows a comparison of the calculated bottom flat wall 
total pressures and the one dimensional analytical solution . 

Both calculations show overshoots and a large undershoot in total 
pressure in the region through the shock. The exit total 
pressures, however, are essentially the same for both the 
calculations and are close to the one- dimensional analytical 
solution. 

The overshoots and undershoots in total pressure arise 
because the pressure used in the momentum equation, an 
interpolated effective pressure, is not used to calculate the 
total pressure. The current code uses the thermodynamic pressure, 
determined from the ideal gas equation of state, to evaluate the 
total pressure. 

Because of the way that the effective pressure is calculated 
in NASA's current f inite- volume code, the shock is smeared out 
over several grid points. One byproduct of this smearing is that 
the maximum Mach number before the shock shown in the 
calculations is lower that the one dimensional value (1.385 and 
1.433 compared with 1.495) and should therefore not be used to 
predict the total pressure loss across the shock. However the 
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total pressure loss through the shock agrees well between the one 

dimensional solution and the calculation. Pt . ./Pt. , . from 

exit inlet 

the one-dimensional solution was .9304 . This compares well with 
the computed value of .931 using the new control volumes and 
reasonably well with the value of .934 using the old control 
volumes. 

The effective pressure is used to stabilize the calculation 

procedure and reduces overshoots and undershoots in static 

pressure and Mach number through the shock. From Fig. 7 , it was 

seen that the local total pressure undershoots considerably 

because of this but that the net total pressure loss through the 

shock is calculated with good accuracy. If the local total 

pressure is calculated from the effective pressure rather than 

the thermodynamic pressure then the total pressure is much better 

behaved as can be seen in Fig. 8. This demonstrates the advantage 

of being consistent by choosing the pressure for use in 

evaluating the total pressure to be the same as the pressure used 

in the momentum equation. Since the effective pressure is equal 

to the thermodynamic pressure at the exit the value of the 

calculated Pt ../Pt. . . is still .931. Fig. 9 shows a 
exit inlet 

comparison of the effective pressure and the thermodynamic 
pressure for this test case . Perhaps the effective pressure 
should be considered the best representation of the "actual" 
static pressure. 

It was attempted to remove the pressure inconsistency brought 
about in the transonic calculations by the way that the effective 
pressure is calculated. The program currently uses Eqs. 2-4 
given in Table 1 to calculate the effective pressure, and for the 
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Table 1 Effective Pressure Calculation 


Pressure Interpolation 

FUP- 1.7 

(1) 

A= (m-i.)*( /l+1 -/ x -i > 

(2) 

/I ■ 

0.0< A. < 0.9 

CFP ] . = (l.-RF)*CFP + RF*(1.-A)*.333*(P x - 

W (3) 

where RF= 0.02 to 0.05 typically and CFP is 

updated 

every 5 iterations 
PEFF = P T , + CFP t 

(4) 

I 1+1 I 

In the limit when convergence is reached, 
PEFF I = P I+1 + (l.-A)*0.333*(P I _ 1 - P I+2 ) 

Relaxation to Ideal Gas 

CFP » (l.-RF)*CFP I + RF*(P I -P ) 

(5) 

PEFF_= P_ . + CFP 

i i+1 I 

(6) 


in the limit when convergence is reached, 
PEFF I = Pj. 
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calculations in Figs. 5-9 the parameter FUP was set equal to 1.7 
(Eq. 1). If Eqs. 5-6 are used instead, as the solution 
approaches a steady state the thermodynamic pressure would equal 
the effective pressure. But this procedure causes large 
overshoots in the static pressure as can be seen in Fig. 10. The 
solution also did not converge and the static pressures shown in 
Fig. 10 are after 2200 iterations. Many different ways of 
applying equations 5-6 were tried for the transonic case but none 
of them got rid of the overshoot problem. However, equations 5 
and 6 could be used to obtain a stable solution if the Mach 
number throughout the duct remained subsonic. 

1.4 THE INFHJENCE OF TRANSVERSE SMOOTHING ON A STEP 

PROFTLE IN A STRAIGHT DUCT 

Transverse smoothing is required in the current Denton 
method with the old control volumes because there are more grid 
points across the duct (unknowns) than there are control volumes 
(equations). Smoothing formulae are used to add non-physical 
"extra equations". Two forms of transverse smoothing are used in 
the Denton code; these are linear smoothing, described in Table 
2, and non-linear smoothing, described in Table 3 and Fig. 11. 
This transverse smoothing of properties causes numerical 
viscosity to be introduced into the solution when large gradients 
in the properties are seen across the duct. 

For the Sajben diffuser' test case, there are no large 
gradients of properties across the duct so you would expect 
little numerical viscosity. This lack of significant numerical 
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Table 2 Linear Smoothing of Flow Properties 


D(J)= (l.-SF)*D(J)+SF*(D(J+l)+D(J-l)) 

2 . 


the variable D at node J is smoothed using this 
equation. The smoothing factor is SF, typically 
0.01 - 0.02. The variables are updated and 
then smoothed. The variables that are 
smoothed are^p , , j> V T ,^rV g , andJe. 


Table 3 Procedure for Non-Linear Smoothing 


1) an average value of a property D is determined from 
the neighboring nodes using linear interpolation. 

AVG(J) (see equation 1 Fig. 11 ) 

2) the difference between the actual and average value of 

a property D at a node is determined and assigned the 
variable name CURVE (J) (see equation 2 Fig. 11 ) 

3) a variable SCURVE is determined from the average of the 
variable CURVE from the neighboring nodes . 

(see equation 3 Fig. 11 ) 

4) the variable D at node I is smoothed using equation 4 

in Fig. 11 . The smoothing factor is SF , typically 0.01-0.02. 

5) this non-linear smoothing procedure results in no 
smoothing added to linearly or parabolically varying 
properties . 
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viscosity is seen by the agreement between the calculations made 
using the new control volumes without smoothing and old control 
volumes with smoothing(see Fig. 6 and 7). As a severe test case, 
calculations were made for a step' profile in inlet properties in 
a straight duct. The geometry can be seen in Fig. 12. A step 
inlet profile of total pressure is specified. The total pressure 
at the centerline is 135 kPa and the total pressure is reduced to 
120 kPa ( pt side /Pt centerline =0 .88 9) at the sides(see Fig. 13) . 
The exit static pressure in the duct is 108 kPa 


(0 .8*Pt centerline ) . 

Fig. 14 shows Mach number profiles at three axial locations 
along the duct for the case where linear smoothing was used (with 
SF=0.02) . The inlet step prof ile{x=0 .Om) is quickly altered 
into a parabolic type prof ile( x=4.0m) . This parabolic profile 
then changes relatively little until the end of the duct 
(x=21.0m). Fig. 15 presents the total pressure distribution 
along the duct. The step profile causes an almost step change in 
the total pressure at the beginning of the duct and then the 
total pressure decreases as in a viscous flow. Fig. 16 compares 
the Mach number profiles at the end of the duct for calculations 
using linear smoothing (SF=0.02) and non-linear smoothing 
(SF=0.02). Non-linear smoothing did not improve the profile. 

Additional calculations were made using the same boundary 
conditions as above but using the new control volumes and no 
smoothing. Fig. 17 compares the inlet Mach number and exit Mach 
number profiles for this case . The improvement over the previous 
results is dramatic. The total pressure distribution has also 
improved especially along the centerline of the duct as can be 
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seen in Fig. 18. These results show conclusively that the 
numerical scheme used to calculate flows with large transverse 
gradients in properties, like those seen in turbulent boundary 
layers, must not have smoothing of properties in the transverse 
direction. 

1.5 SHOCK LOSSES IN AN INCLINED DIFFUSER 

One of the possible sources of inaccuracy in the calculation 
of total pressure in transonic compressor calculations is due to 
the possibility of smoothing through the shock due to a inclined 
flow path. In the Sajben diffuser calculations for Figs. 6-10 
the property gradients across the duct were not large and 
transverse smoothing did not add noticeable numerical viscosity 
into the calculations. However, if the diffuser were inclined at 
an angle, the shock would become oblique to the grid and would 
introduce large transverse gradients in properties there. Fig. 

19 shows Sajben's diffuser inclined at an angle of 40 degrees 
with respect to the horizontal axis. This geometry would cause a 
normal shock to cross about 4 transverse grid lines. The Sajben 
geometry used previously( see Fig. 5) has been extended with 
constant area sections added on to the inlet and exit of the 
duct. The same inlet and exit boundary conditions were specified 
as in the previous Sajben diffuser calculations. The input file 
used in these calculations is in Appendix A. 

The effective pressures along the flat wall are plotted in 
Fig. 20 for the old and new control volumes, with and without 
smoothing respectively. The wall pressures using the new control 
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volumes were determined by linearly extrapolating from the 

interior points. The effective pressures for the old control 

volumes are displaced downstream of the effective pressures using 

the new control volumes. Fig. 21 presents the total pressures 

along the bottom wall calculated using the effective pressures. 

The exit total pressures are approximately the same(0.941 for new 

C.V. and 0.944 for old C.V.) for both calculations but the local 

total pressures from the calculation using the new control 

volumes are better behaved. In both cases, however, the exit 

total pressure ratios of 0.941 and 0. 944 are di ff e rent from the 

exit total pressure of 125.7 kPa (Pt . ,/Pt. . .=0.931) 

exit inlet 

calculated for the horizontal sajben diffuser. 

Further comparisons of the inclined and horizontal results, 

using the old and new control volumes, are shown in Figs. 

22 ,23 ,2 4, and 2 5. The minimum static pressure for the inclined 
calculations is greater than that obtained from the horizontal 
calculations. The total pressure behavior is also poorer for the 
inclined calculations when compared with the horizontal solution. 
The total pressure losses through the shock for the inclined 
calculations are approximately 15% less than those for the 
horizontal calculations. The oscillations in total pressure at 
the exit of the inclined diffuser (with the old control volumes) 
are perhaps the result of the long thin control volumes that are 
seen there. 

It was difficult to obtain a converged solution for the 
inclined geometry with either control volume. So the results 
presented here are after 1000 iterations. At this point in the 
calculations the maximum error in mass flow rate with the old 
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control volumes is 0.2% and the maximum change in axial velocity 
was .00009 of the mean flow velocity. For the new control 
volumes, the corresponding error is 0.2% and the corresponding 
change is .00019. 

1.6 EFFECTIVE DENSITY CALCULATIONS 

The method presented in reference 1 uses a different 
procedure to update the pressure and density. The pressure is 
updated directly from the continuity error and then the density 
is updated using the ideal gas equation of state. This procedure 
was necessary because of the multi- volume approach used in the 
boundary layer region. This updating procedure was also 
implemented in the three dimensional version of the finite volume 
code at NASA Lewis. 

If the density is updated in the calculations such that the 
effective density becomes the actual density at a node the same 
overshoot phenomena in static pressure and Mach number appears 
here as did when the effective pressure was used (see Fig. 10). 
The solution would also not converge. Therefore an effective 
density which does not use the actual pressure at a node, but 
uses an interpolated pressure, is used and is described in Table 
4. This effective density reduces the overshoot problem and 
results in a stable solution. The following results use this 
effective density. Just as the interpolation procedure before 
introduced an inconsistency between the thermodynamic and 
effective pressures, the ideal gas equation of state is not 
satisfied completely when the effective density is used. 
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Table 4 Effective Density Calculation 


CFPj = (l.-RF)*CFP + RF/3.*(P I+1 - P ) (1) 




= (Pj + CFP I ) 
L I+1 


( 2 ) 


R*T, 


in the limit when convergence is reached 


1+1 P I + P I+l~ P I-2 


RT 


1+1 
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Fig. 2 6 compares the bottom flat wall static pressures for 
calculations of transonic flow in Sa.jben's diffuser (Fig. 5) 
using the effective pressure method and the effective density 
method with the same boundary conditions as in our previous test 
cases and the new control volumes (the effective pressure is 
shown for the effective pressure method) . Fig. 27 compares the 
total pressures for these two cases. The effective density 
method gives a much more uniform total pressure upstream and 
downstream of the shock; there are no overshoots in total 
pressure when the effective density method is used. 

To obtain a stable solution using the effective 
density, it was found necessary to assume a constant total 
temperature rather than use the energy equation in its full 
form. It may be that an interaction between the continuity 
error and the energy equation was responsible for this 
instability. 

1.7 CASCADE GEOMETRIES 

Some additional calculations were made, using the new control 
volumes and the effective pressure method, on some simple cascade 
geometries. The purpose of these runs was to check out the 
periodicity condition and the treatment of the leading and 
trailing edges which were discussed earlier. The geometries are 
shown in Figs. 28 and 29. The inlet total pressure was 101.352 
kPa, the inlet total temperature was 288.166 K, and the exit 
static pressure was 85.44 kPa. Fig. 30 shows the Mach numbers 
calculated along grid lines which are closest to the pressure and 
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suction surfaces. A copy of the input data is in Appendix A. 

The only difficulty encountered in these cascade flows was a 
problem with the flow properties oscillating at the inlet about 
some mean value between adjacent grid points. This was due to 
the periodicity condition and the absence of smoothing to damp 
out these oscillations. This problem was found to occur only 
when there were an even number grid points across the inlet. An 
odd number of transverse grid points seems to decouple the 
odd-even oscillations. 

1.8 CALCULATIONS FDR AN INLET GUIDE VANE 

Fig. 31 shows the geometry and grid for an inlet guide vane 
that was my final test for the summer work. The total pressures 
along the streamline closest to the suction surface are presented 
in Fig. 32. The new control volumes are used for both 
calculations, one uses the effective pressure method and one uses 
the effective density method. The total pressure distributions 
in Fig. 32 , which use the new control volumes ,can be compared 
with the results ,from the Denton code using the old control 
volumes, shown in Fig. 33. Both the calculations using the 
effective pressure method have large oscillations in total 
pressure around the leading edge whereas the calculations using 
effective density method give better total pressure behavior. 

The effective density method is much better at calculating the 
total pressure than the effective pressure method even when the 
effective pressure is used to calculate the total pressure as was 
done for the results shown in Fig. 32 . 
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Pigs. 3 4 and 3 5 present Mach numbers along grid lines 
through the inlet guide vane for calculations using the effective 
density method with the new control volumes and the effective 
pressure method with the old control volumes /respectively . It 
was found necessary to use a small amount of smoothing ( SF=0 .00 5) , 
even using the new control volumes ,for this inlet guide vane. 
This need for smoothing is perhaps the result of an incorrect 
treatment of the pressure at the trailing edge. 

SUMMARY 

The work at NASA Lewis Research Center . involved using and 
modifying the Denton f inite- volume time-marching code. The 
results can be summarized as follows. 

1. The new control volumes developed at V.P.I. & S.U. can 
be extended to a three dimensional geometry. For duct 
flows the new control volumes require no transverse 
smoothing. 

2. The results for Sajben's diffuser were essentially the 
same for both control volumes. 

3. If the effective pressure is used to calculate the total 
pressure the total pressures are much better behaved 
through the shock. 

4. Even though the total pressure is incorrectly calculated 
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in the shock, the overall total pressure loss through the 
duct is calculated accurately. 


A. One dimensional solution 


B. Old control volumes 


C. New control volumes 


Pt exit //pt inlet 


0.93 0 4 


0 . 93 4 
0.931 


5. An interpolated effective pressure is needed to stabilize 
the solution. 

6. For calculations where large transverse gradients in 
properties are observed, transverse smoothing cannot be 
used if an accurate solution is to be expected. 

7. The new control volumes with no transverse smoothing 
allow calculations with large transverse gradients in 
properties without significant numerical total pressure 

. losses. 

8. Good convergence was not obtained for the inclined Sajben 
with either control volume. 

9. Calculations which use the effective density method, 
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developed at V.P.I. & S.U. , rather than an effective 
pressure method give a much .more uniform total pressure 
upstream and downstream of the shock. Overshoots and 
undershoots were a factor of 10 smaller with the 
effective density method. 
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Fig. 2 New Control Volumes in Non-Uniform Grid 



Fig. 3 Periodicity Condition 
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Fig. 5 Geometry and Grid for Sajben's Diffuser Calculations 
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Fig. 11 Non-Linear Smoothing 
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Sajben's Diffuser Inclined at 40° Angle 
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HORIZONTAL AND INCLINED 
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FOR SAJBEN DIFFUSER-NEW CONTROL VOLUMES 
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Input File for Inclined Sajben's Diffuser 
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NOMENCLATURE FOR PART 2 


A area 

a Q ,a 1 ,a 2 pressure interpolation coefficients, Eq. 51 

c speed of sound 

e =0 for pressure update method 

=1 for density update method 
i grid index 

M Mach number 

A mass flow rate 

p static pressure 

R gas constant 

T static temperature 

u, n velocity in x direction, velocity vector 

Vol volume 

x 1-d coordinate 

&t time step 

St time step for continuity 

p density 

y ratio of specific heats, c /c 

p v 
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PART .2 

1-D STABILITY ANALYSIS OF DENSITY-PRESSURE RELATIONS 
USED IN THE COMPUTATION OF TRANSONIC FLOW 

2 . 1 Background 

When calculations of l-d or 2-d choked flow with a shock were 

* 

attempted with equations that were relaxed to perfect gas (i.e. 
so that when converged, the ideal gas equation of state would be 
satisfied with the same density and pressure as used in the 
momentum and continuity equations) convergence was not obtained. 
Eventually, as more and more iterations were taken wobbles 
appeared in the pressure solution which grew and continuity 
errors grew worse instead of better. 

The following analysis explains the cause of the instability. 
Further analysis then shows the stability of the 3 point 
interpolation scheme used for the calculation of effective 
pressure. Still further analysis suggests a Mach number dependent 
interpolation scheme. 

2 . 2 l-D Flow Example 


~r — 
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1 
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We are seeking a l-d steady flow solution for continuity and 
momentum 


Vpu = o 

(1) 

Vpu li = - Vp 

(2) 


for a perfect gas with constant total temperature. 


* Using either an effective pressure (Denton) or an effective 
density (Nicholson/Moore) finite-volume time marching method. 
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2 . 3 Continuity 


Continuity between grid points i and i+1, in discretized 
form, when a converged solution is obtained, is simply 


f f f f 

> . , i u . , , A. , , - p . u . A. = 

l+l l+l l+l i ii 


( 3 ) 


where superscript f stands for final. 

Now consider an intermediate solution, p and u, which does not 
satisfy continuity and changes, 6p and 8u, so that continuity is 
satisfied. Then 

(p i+l +8p i + l ,(U i + l +8u i+l )A i-fl " <P 1 + 6P i ><u.+ Su i> A i = ° (4) 

Rearranging, 

P • . i A. 8u . - p.A.Su. + U ■ , . Ai , 8p . , - u.A.Sp. 

l+l l+l i+l ii l l+l l+l l+l ill 


= p i u i A i ' p i+l u i+l A i+l + 8 p i 8u i A i ‘ 8p i + l 8u i + l A i« 


( 5 ) 


The first two terms on the right hand side represent the current 

2 

continuity error and the last two are of order 8 and so will be 
negligible when the computation is nearly converged and 8p<<p and 
8u< <u . Therefore we may write this equation as 

p i-*-l A i-t-l 8u i+l - p i A i 8u i + u i + l A i-n 8p i + l - u i A i 8p i 


= ft . + small 

error , i 


( 6 ) 


In the density update time marching calculation procedure 
(Denton), when there is a continuity error, the density on the 
downstream side of the control volume is changed. 


8 p 


i+l 


wO * 

e r ro r 


^t/Vol.. 


( 7 ) 


The density change affects continuity directly, but it also acts 
through the perfect gas equation to change the pressure which acts 
through the momentum equation to change the velocity. 


In the pressure update time marching calculation procedure 
(Nicholson/Moore), when there is a continuity error, the pressure 
on the upstream side of the control volume is changed. 



^error 


. 8 t RT/Vol . 

1C l 


(8) 
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This pressure change acts through the momentum equation to change 
the velocity and through the perfect gas equation to change the 
density. 


2.4 Momentum 

The .steady state momentum equation discretized over the 
control volume between points i and i+l is 

(PUA , 1+1 U 1+1 - (puAJjUj = -fP i+1 A ltl - - PgCA^j-A^J (9) 

where p is the pressure acting on the sides of the control 

D 

volume. Traditionally, 

P s ’ ( Pi.l + Pi >/J < 10) 


therefore. 


(puA). +1 u. +1 - (puA).u. = -(p i+1 - pj.) (A i+1 +A.)/2 (11) 

We may write this as 

( 12 ) 




*'i u i " * ( Pj+i ~ Pilvolj/Sx. 


where lti=puA is the mass flow rate, Vol i =6x i (A i+1 +A i ) /2 is the 
volume of the control volume and 8x i =x i+i _x i is the grid spacing. 

Eq. 12 may be rewritten as 

<u i+i _u i )( *i+i + *i ,/J " < *i+rV <u i+i +u i )/2 ' - ( Pi tl -Pi' Vol i /8,! i 

(13) 

or 

" (p i+r p i )Vol i /8x i (14) 


»(u itl - u i) + u !» error>i 


In the Nicholson/Moore method the continuity error term is 
omitted and the change in velocity on the downstream side of the 
control volume is proportional to the momentum error. 


& u. +1 = t - <P i+1 +8p i+1 -p i -5p i )Vol i /5x i 

- ft(u i+1 -u i ) ] 5t/ (p i+1 Vol i ) (15) 

where 6p is the change in pressure calculated from the continuity 
error. In the Denton method, the continuity error is not omitted 
in the momentum equation and the change in pu is calculated from 
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the momentum error. 


8(pu). +1 = u i+1 Sp. +1 + P i+1 8u. +1 

= l -<P i+1 + &P i+1 -P i -&P i )Vol i /&x i 

- “ terror, i 1 5t/Vol i (16) 

Taking the mean velocity u approximately equal to the velocity on 
the downstream side u^ + 1 , we may subtract u times Eq. 7 from Eq. 

16 and so for the Denton method we get 


Su 1+1 - [ - <P 1+1 +»P i+1 -P 1 -»P i >Vol 1 /8x i 
- ft(u i + 1 -u i ) 3 St / (p i+1 Vol i > 

the same as for the Nicholson/Moore method. 


(17) 


If at the beginning of a time step the momentum equation is 
balanced except for the continuity error. 


Mu i+1 -u L ) = -<P i+1 -P i >Vol./8x i 
then for both methods we have 

Su i+1 = (8 Pi - S Pi + l> st/<p i+1 8x.) 


(18) 


(19) 


In general , in the density update (Denton) method, the time 
step is calculated from the CFL condition 


8 t = 8x/(u+c) 


( 20 ) 


where c is the speed of sound. In the pressure update method the 
time step for momentum is obtained from the coefficient of u i+1 

in the steady flow equation so that 


8 t = 8x/u. 

We may combine these two equations by saying 
8t = 8x/(u+ec) 


( 21 ) 


( 22 ) 


where e=l for the density update method and e=0 for the pressure 
update method. Combining Eqs . 18 and 15 then gives 

p i+l 8u i+l = (5p i ' 5 Pi + l )/(u+ec) i + i* (23) 


2.5 



2 . 5 Change in Continuity for One Time Step 

The left hand side of Eq. 6 may be used to evaluate the change 
in continuity for one time step. Substituting Eq. 23 into this 
expression to eliminate p8u yields 

A i+l (8p i _8p i+l )/(u+ec) i+l " A i (8p i-l~ 8p i )/(u+ec) i 

+ U . . , A. . 8p . , - u.A.Sp. = ft . .. (24) 

l+l i+I * 1+1 ill change , i 

Rearranging to order the coefficients of the 8p's and 8p's 


- A^/(u+ecK 

6p i-l 


+ [A^ +1 / (u+ec) i+j+A^/ (u+ec) 

8p i 

- U i A i 8p^ 

- A 1+1 /(u+ec) 1+1 

6p i + l 

+ “i.iVi 6p i.i 


= ft , . (25) 

change, l 

For stability we require that the change in continuity be of 
the same sign as the continuity error. Note that this is a 
necessary condition for stability but may not be a sufficient 
condition to ensure stability. 


2 • 6 Stability of Density Update Method Using Perfect Gas 


For an intermediate solution where there is a continuity error 
only between i and i+1, Eq. 7 yields 


8p i+l = ft error,i 8t/Vol i 


8p i 


8 P 


i-1 


= 0 
= 0 


(26a) 

(26b) 

(26c) 
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and from perfect gas, assuming temperature changes over the time 

* 

step are negligible , 


5p. . = 8p.^,RT = ft .RT 5t/Vol. 

*i+l 1+1 error, i 1 

5p . = Sp.RT =0 

*1 1 

8p i-l = 5p i-l RT = °* 


(27a) 

(27b) 

(27c) 


For this case then Eq. 25 reduces to 
r - a i+ 1 RT/(u + c ) 1+1 * a i+1 » i+1 ) » etror ,i 8fc/Vol i ’ * 


(28) 


Since for stability we require ft and ft , 

J ^ error change 

sign, we must have 

[ - RT/ (u+c) . , + u.^, ] A. , > 0. 

l+l l+l l+l 


change, i 
to have the same 

(29) 


Substituting c fy for RT yields 


or 


- c / Ir (u+c) i+1 ) + u i+1 > o 


r u. +1 (u+c) itl > c 2 . 


Evaluating for r = 1.4 we need 
u > o.48c. 

Thus for low Mach number flow this density update method is 
unstable. 


(30) 

(31) 

(32) 


2 • 7 Stability of Pressure Update Method Using Perfect Gas 

For a continuity error only between i and i+1, Eq. 8 yields 


6 Pi " *error,i 5t c RT/Vo:1 i (33a > 

6p » 0 (33b) 


* If at this point the alternative assumption was made that the 
changes were isentropic, then 8p = yRTSp and less conservative 
stability criteria would be obtained, equivalent to setting y=l.o 
in the following analysis. 
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0 


(33c) 



and from perfect gas assuming temperature changes over the step 
are negligible 


Sp. = SPj/RT 

6 » i+ i • 5 p 1 *i /r ' 1 - 

Therefore, in this case, Eq. 21 becomes 


(34a) 

(34b) 


(A. . , /u. , +A./U. -u.A./RT)A .St RT/Vol . = ft . .. (35) 

i+l l+l li li error, l c l change, l 

For ft error and to have the same sign, we require 

A. .... / U • . , + A./u. - u.A./RT > 0. (36) 

1+1 i+l 11 11 

Assuming the values at i are approximately the same as the values 

2 

at i+l, and again using c ty for RT yields 


2/u - 

yu/c 2 > 0 

(37) 

or 



2c 2 > 

ru 2 . 

(38) 

For r = 1.4 then we 

need 



u < 1.2 c; (39) 

thus for high Mach numbers this pressure update method is 
unstable . 


2 . 8 A Downwind Effective Pressure or Upwind Effective Density 

If an inconsistency in the pressure-density relation is 
introduced such that the pressure used in the momentum equation is 
offset by 1 grid point from the density used in the continuity 
equation, the equation of state may be written as 

Pi - P i+ iRT < 40 > 

In a density update method this may be viewed as an effective 
pressure evaluated downwind of its point of use in the momentum 
equation. Similarly, in a pressure update method the equation 
represents an effective density evaluated upstream of its point of 
use in the continuity equation. 
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For both the density and pressure update methods, Eq. 21 now 
becomes 


- (A. / (u+ec) . + U.A. /RT) (StRT/Vol. . ) 
1 111 1-1 


ft 


error, i-l 


+ (A. . .. / ( u+ec) . , , +A. / (u+ec) . +u. , A. . /RT) (StRT/Vol. ) ft 

l+l l+l l l l+l l+l l error, l 


- (A^ +1 / (u+ec) i+1 ) ( StRT/Vol i+1 ) 


ft 


- ® h 


ange, l 


error , i+l 
(41) 


From this equation we can see that the coefficient of ft .. .is 
^ er ror , l 

always positive and so this pressure-density relation passes the 

simple stability criterion (ft , . same sign as ft . ) for 

r J change , i error , i 

all Mach numbers. Note also that the coefficients of ft . , 

error, i-l 

and ft . .. are of opposite sign to the coefficient of ft . 

error, i+i error , i 

and in general of smaller magnitude; this further assures the 

stability of Eq. 40. 

While the pressure-density relation, Eq. 40, is stable, 
testing has shown that it results in poor shock capturing as the 
calculated shock is spread over numerous grid points. Fig. l 
shows the calculated and theoretical pressure distribution for a 
1-d calculation with a nominal shock Mach number of 1.45. 
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2.9 Stability of 3-p oint Interpolation for Effective. Pressure 

One of the pressure-density relations used in the density 
update method is a three point interpolation of density to obtain 
the effective pressure. For approximately uniform T we may write 
this as 

Pi = (p i + l " (1/3)(p i+ 2 - p i-l> > RT (42) 

so that for the change in p we have 

8p i = (8 Pi+i " (1/3)6 Pi+2 + (1/3 )5p i _ 1 )RT. (43) 

Substituting into Eq. 25 and neglecting variations of A, u 
and c with i, we obtain 

- (Ac 2 / (y(u+c) ) ( 8 Pi -<l/3)6p i + 1 + ( 1 /3 ) S p ± 2 ) 

+ (2Ac 2 / (r (u+c) ) ( 8p -<l/3)8p. , +(l/3)8p. ) - uA 8p. 

- (Ac 2 /(t<U+c) ) ( Sp. +2 -(l/3)8p. +3 +(l/3)5p.) + uA Sp i+1 


= flL 


(44) 


changed 

Collecting terms and substituting the Mach number, M for u/c 

+ Ac/ ( 3 y ( M+l ) ) 

- 5 Ac/ (3 y (M+l ) ) 

+ Ac ( 7 / (3 y (M+l ) ) +M) 

- Ac ( 4 / ( 3 y ( M+ l ) ) +M ) 

+ 2Ac/ (3 y (M+ l ) ) 

- Ac/ (3 y (M+ l ) ) 8 p . „ = ffl 


5p i + 3 
5p i + 2 
6p i + l 
8p i 
6 p i-l 

8 p . _ 
* 1-2 


change, i 


(45) 


we can see that the change in continuity for the control volume 
between i and i+l is now dependent on the change in density at 6 
grid points. 

Eq. 45 passes the first simple test for stability; the 
coefficient of 8 P^ +1 is positive for all Mach numbers and so a 

continuity error, rtt ^ . , will result in a change in 

error, l 

continuity, lli^ of the same sign. However, since the 
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coefficients of &p^ +3 and Sp^,^ are also positive, it is 

appropriate to apply a more sophisticated stability criterion. The 
criterion we will apply is: 

the center point coefficient must be greater than the sum of 
the other positive coefficients. 

Coef _ > Sum Coef^ (46) 

center + 

Applying this to Eq. 45, we require 

Ac (7/ ( 3y ( M+l ) ) +M) > 3 Ac/ (3 y (M+l ) ) (47) 

which is always true. Thus, Eq. 42 should be stable for all Mach 
numbers. The experience of Denton and other users of his code 
confirms this. 


2.10 stability of 3-point Interpolation for Effective Density 

A similar analysis can be done for the pressure-update 
effective-density method using a three point interpolation of 
pressure to obtain the effective density. 

P i + 1 - t Pl + (l/3Hp i(1 -p._ 2 ))/M (48) 


Substituting the change in density 


8p i + 1 = Up. + (l/3>5p. +1 


(l/3)6p i _ 2 )/RT (49) 


into Eq. 25, gives 


+ 

(A/c) (-1/M + 

yM/3 ) 

8 p i+ i 

+ 

(A/c) ( 2/M + 

2yM/3 ) 

8 Pi 

- 

(A/C) ( 1/M + 

yM ) 

S Pi-l 

- 

(A/c) ( 

yM/3) 

6Pi- 2 

+ 

(A/C) ( 

yM/ 3 ) 

S Pi-3 


.a. 

m change, i 


(50) 


The center point coefficient is the coefficient of Sp^, siace 

this is proportional to * error ^ in the pressure update method. 

In Eq. 50, the coefficient of &p^ is positive and greater than 

the sum of the other positive coefficients; therefore Eq. 48 
should be stable for all Mach numbers. 
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2.11 Mach Numb er Depe ndent, Interpolation Formula 
for Effective Density 

In section 2.7 it was indicated that when the Mach number is 
low, the pressure update method is stable with the ideal gas 
equation of state satisfied at each grid point. Since this is the 
correct pressure-density relation for ideal gases it should be 
used where feasible. In this section we will start with a 
generalized pressure interpolation equation for effective density 

p i + i = l ?i + a o (p i + r p i ) + a i (p i+ r p i-i )/2 

+ a 2 ( Pi+l" P i-2 ) /3 3/RT (51) 


and seek Mach number limitations to a Q , aj and a 2 using criterion 

(Eq.) 46. Comparing equations 40 and 51, the upwind effective 
density corresponds to a Q =a 1 =a 2 =0, the 3 point interpolation, Eq. 

48, corresponds to a o =a 1 = 0 » a 2 = l, and ideal gas to a o =l, a 1 =a 2 =0. 


Substituting 


8p 


i + i 


[(1_a o )8p i + < a Q + a 1 /2 +a 2 / 3 > 8 P i+1 
- - (a 2 /3)8p i-2 1/RT 


(52) 


into Eq. 25 and rearranging in terms of the ceofficients of each 
8p, a Q , aj and a 2 , yields 


(A/c) ( (-1/M + rMa Q + 

+( 2/M +yM - 2yMa o - 
+(-l/M -yM + YMa Q - 
+ ( + 
+ ( 

_ JX 

change, i 


<YM/2)a 1 
( yM/ 2) 

( yM/ 2 ) a^ 
(YM/2)a x 


+ (yM/ 3 )a 2 > 8 P i+1 

- ( y M/ 3 ) a 2 ) 

> spj.j 

- (YM/3)a 2 ) 8 P ^_2 
+ (YM/3)a 2 > SPi.3 

(53) 


Let us first consider the case when a 1 =a 2 =0 and find limiting 

values of a . From Eq. 51, it is obvious that we should consider 
o 

only values in the range 


0 < a < 1. 
— o 


(54) 
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( 55 ) 


The coefficient of 8p^ is positive when 

2/M + yM - 2 Y Ma 0 > 0. 

This gives a limit on a which is a function of Mach number, 

3 o 

a < l/(yM 2 ) + 1/2. (56) 

o 

But the coefficient of s P i+1 is positive when 

-l/M + yMa. > 0, or m 2 > l/(ya ) (57) 

o o 

In this region, from Eq. 46, we require 

2/M + yM - 2yMa o > -1/M + yMa Q (58) 

or 

a 0 < l/(yM 2 ) + 1/3 . (59) 

Valid values of a„ based on these criteria are shown as a function 

o 

of Mach number in Fig. 2. 

Let us next consider limiting values of a x when a Q and a 2 are 
zero. From Eq. 53 the coefficient of Sp^ 

2/M + yM - yMaj/2 > 0 (60) 

is positive for all Mach numbers in the range 

0 1 * 1 1 1 . (61) 

The coefficient of 8 P^_ 2 is positive for all M and the 
coefficient of 8 Pj_ + 1 is positive when 

-1/M + yMa t /2 >0 or M 2 > 2/(ya 1 ). (62) 

2 

For M < 2/(ya 1 > we then require the coefficient of 5p^ to be 
greater than the coefficient of 8 P^_ 2 » 

2/M + yM - yMa x /2 > yMaj/2. (63) 
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. 2 
With a 2 < 1, this is always satisfied. For M > 2/(Ya x > we 

require the coefficient of Sp^ to be greater then the sum of the 

coefficients of 8 P^_ 2 an<3 8 P^ +1 » 

2/m + rM - rMaj/2 > yMaj -l/M (64) 

or 

a 1 < 2/(yM 2 ) + 2/3. (65) 

2 

Thus we can have a 1 =l up to M = 6/y or up to M=2.07 for y=l.4. 

Fig. 3 shows the valid range of a^ based on these criteria. 

We now consider combinations of a Q , a 2 and a 2 . In particular 
if 

a Q + a x + a 2 = 1. (66) 

the interpolation scheme is second order accurate. (See Appendix 

A.) 

For Mach numbers less than 2, a 1 =i is stable. Therefore, for 
M < 2, we will choose 

a 2 = ° 

a Q + a x = l. (67) 

From similar stability analyses to those already given 

a Q < 2/(yM 2 ) - 1/3 (68) 

should be stable for M <. 2 . 

For M > 2, we will choose 

a 0 = 0 

a 2 + a 2 = 1. (69) 

The stability analysis suggests acceptable values of a 1 are 

a 2 < 0.4 + 3.6/(yM 2 ) . (70) 

The stability criteria, Eqs. 68 and 70, are shown on Fig. 4. 
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A set of equations for a Q , a 2 and a 2 , which satisfy Eq. 66 and 

so give second order accurate interpolation, and which satisfy 
Eqs. 6 7-70 so that they satisfy the stability criteria, have been 
selected. These are: 


for M < 2 


for M > 2 


a Q = (0.8/3) (4/M 2 -1) 



(71) 


(72) 


These Mach number dependent formulations for a Q , a x and a 2 are 

shown in Fig 5. These equations are tested in Part 3 where they 
are referred to as tke M&M formula. 
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STABILITY LIMITS. A1=A2=0 



0.00 0.20 0.00 0.60 0.80 >.00 1.20 I.WO 1.60 1.80 2.00 2.20 2. VO 2.60 2.80 3.00 

MACH 


Fig. 2 Acceptable values of a^ as a function of Mach number 
based on Eqs. 54, 56, 57, and 59 (for y = 1.4). 


2.17 



1 . OQ 


STABILITY LIMITS, A0=A2=0 



0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 

MACH 


Fig. 3 Acceptable values of a^ as a function of Mach number 
based on Eqs. 61, 62, 63, and 65 (for y = 1.4). 
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STABILITY LIMITS, flO+Rl = 1 . A2=0; R0=0. R1 +R2= 1 



0.00 0.20 O.UO 0.60 0.80 1.00 1.20 l.MO 1.60 1.80 2.00 2.20 2.10 2.60 2.80 3.00 

MACH 


Fig. 4 Stability limits for a Q when a^ = 0 and a Q + a 1 = 1, 
and for a^ when a^ = 0 and a^ + a^ = 1. 
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MRCH NUMBER DEPENDENT fl’S WITH flO + Rl+R2=l 


O 



0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 2.80 2.00 2.20 2.40 2.60 2.80 3.00 

MflCH 


Fig. 5 


M & M Mach number dependent values (Eqs. 71 and 72) 
for the coefficients in Eq. 51. 
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PART. 2 


APPENDIX A. TRUNCATION ERROR OF PRESSURE INTERPOLATION EQUATION 

The truncation error of the interpolated pressure used to 
calculate the density in Eq. 51 may be determined using Taylor 

series analysis. The interpolated pressure p e is given by 


i+l 


■ Pi + a o'Pi + 


1-Pi* 


a l ( Pi*l-Pi-l> /2 + VPi.l-Pi-2 


) / 3 


and to determine the accuracy of p e we will look at the magnitude 
of p e ^ +1 -p^ +1 . With grid spacing h, arid expanding about i+l, we 
have 


Pi-2 

= P “ 

3hp ' + 9(h 2 /2)p" 

- 0(h 3 ) 

Pi-1 

■ P - 

2hp ' + 4(h 2 /2)p' ' 

- 0(h 3 ) 

Pi 

- P “ 

hp' + (h 2 /2)p' * 

- 0(h 3 ) 

Pi + 1 

= P 




Therefore, 

p e ^ +1 ~p^ +1 = hCaQ+aj+aj - ! )p ' “ (h 2 /2) ( a Q +2a 1 +3a 2 -l )p ' ' + 0(h 3 ). 
And if 


+ a. + a. 


1 , 


2 

then the difference between p^ and p is of the order of h , so 
e e 

that p is a second order accurate approximation for p. 
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a o' a l' a 2 

k = c p /c v 
M 
p 
R 
T 
x 

6 X 
P 

Superscript 

e pressure used to evaluate the effective density 

Subscript 
t total 


pressure interpolation coefficients, Eq. 2 

grid index in flow direction 
ratio of specific heat capacities 

Mach number 
pressure 
gas constant 
static temperature 
axial distance, Eq. 1 
grid spacing 
density 
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PART 3 

1-D COMPUTATIONAL TESTS OF SHOCK CAPTURING USING PRESSURE 
INTERPOLATION FORMULAE TO CALCULATE EFFECTIVE DENSITY 

3 . 1 Denton's 1-D Nozzle for Testing Shock Capturing 

Denton [1] has tested shock capturing with his finite-volume method 
in a convergent-divergent nozzle (see Fig. 1) designed to produce a 
linear variation of Mach number with distance for 1-D isentropic flow. 
The equation for the Mach number variation with distance is 

x = 10. + 45. (M - 1) (1) 

Denton considered flow between x = 1, M = 0.8 and x = 46, M = 1.8; the 
throat, M = 1.0, is at x = 10. He used three back pressures with 

p exit^ p t inlet = 0.80, and 0.75, respectively. The theoretical 1- 

D solutions for these three flows are shown in Fig. 2. The maximum Mach 
numbers, just upstream of the shock, are 1.267, 1.455, and 1.578, re- 
spectively; this is a range of shock Mach number typical of turboma- 
chinery flows. 

We have used these three pressure ratios for Denton's 1-D nozzle to 
test shock capturing with three of the pressure interpolation methods 
discussed in Part 2. 


3.2 Effective Density Method 

This annual report (Parts 1, 2, and 3) includes the results of 

calculations made with three different methods: 

(a) The Effective Pressure Method as currently programmed in Denton's 
code at NASA Lewis. 

(b) The Effective Density Method incorporated into Denton's code at 
NASA Lewis by S. Nicholson; this method uses the same time steps 
for the continuity and momentum equations. 

(c) The Effective Density Method developed at VPI&SU with different 
time steps for continuity and momentum [2,3]; this is incorporated 
in the Nicholson/Moore time-marching codes at VPI&SU. 


These methods are outlined to show their similarities and differences in 
Table 1 . 


In Part 1, results from methods a and b were presented and com- 
pared. The stability analysis in Part 2 was applied to methods a and 
c. Here in Part 3 the calculations are performed using method c. 


3.3 Pressure Interpolation Schemes 


The effective density methods (b and c in Section 3.2) use an 
interpolated approximation for the pressure in the evaluation of the 
density. A general form of the interpolation formula considered in this 
report is 


P i+1 P i + a o^ P i+l P i ) + 2 ^ P i+1 P i-1^ + 3 (P i+l ~ P i-2 ) ^ 
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Table 1. Outline of Effective Pressure and Effective Density Calculation Methods. 
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MOMENTUM (IN VISCID) 
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constant 



( 3 ) 


and this Is used to evaluate the density as 

P e 

_ 1+1 

P i +1 “ RT 1 + 1 ' 


This general form is a linear combination of a single-point interpola- 
tion, P 1+1 - P i} a 2-point interpolation, P i+1 - Pj_j , and a 3-point 
interpolation, P^ + j - Pi-2* The single-point interpolation, of course, 
really gives the correct perfect gas equation and involves no approxi- 
mation. 


The coefficients a Q , a^ , and a 2 are here taken to be constants or 
functions of Mach number. Combinations, including individual terms or 
pairs of terms, for which the sum of the coefficients 

a Q + a : + a 2 = 1 

are second order accurate, as shown in Appendix A, Part 2. 

Correct Perfect Gas Equation (a Q = 1, a^ = 0, - 0) 

This scheme has the advantage that it involves no interpolation or 
approximation for the pressure. Experience has shown (see Part 1) that 
it is stable for subsonic flow. But the stability analysis of Part 2 
shows that for Mach numbers above about 1.2 this scheme becomes un- 
stable. Thus it could not be used for the test cases of Denton's 1-D 
nozzle. 

These observations about the use of the correct perfect gas equa- 
tion are in agreement with Denton’s findings for his scheme B [1], In 
that method changes of density were sent to the upstream corners of the 
element, which is equivalent to our sending pressure changes upstream. 
The method "proved stable, without any correction factors or damping, at 
low Mach numbers but instability was found to develop at Mach numbers 
around unity and above." 

2- Point Interpolation (a Q = 0 , = 1, a 2 = 0) 

The stability analysis of Part 2 shows the 2-point scheme to be 
stable for Mach numbers up to about 2.0. Use of a 2-point scheme or a 

3- point scheme has been suggested by Denton in his recent ASME and AGARD 
Lecture Notes [4,5]. 

3-Polnt Interpolation (a Q = 0, a^ = 0, a 2 = 1) 

3-point schemes have been shown in Part 2 to be the most conserva- 
tive (in terms of stability) of the schemes considered in this report. 
Perhaps for this reason, such a method is used to stabilize the current 
NASA version of the Denton code. Both 3-point and 2-point schemes 
provide second order accuracy for a continuously changing pressure; they 
give correct interpolated values for linear variations in pressure 
(assuming equally spaced grid points). 
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M&M Mach Number Dependent Interpolation 


The advantages of the three schemes just considered are: 

(1) the accuracy and stability of the perfect gas equation for subsonic 
flow; 

(2) the stability of the 3-point interpolation at Mach numbers greater 
than 2.0; 

(3) the stability and reduced smearing of properties of the 2-point 
interpolation at supersonic Mach numbers up to 2.0. 


These advantages have been combined in a single Mach number depen- 
dent interpolation scheme in Part 2. In this method 
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The values of the three coefficients are shown graphically in Fig. 3. 
Note, once again, that the sum of the coefficients is equal to one for 
all Mach numbers, so that this scheme is also second order accurate. 
For the calculations presented here M was taken as the larger of M on 
the upstream side or downstream side of the control volume. 


3.4 Computational Tests of Three Pressure Interpolation Schemes 

Of the four schemes just considered, three are stable in the Mach 
number range 1.0 to 2.0. These are the 2-point, 3-point, and M&M inter- 
polation methods. In this section, results of shock capturing with 
these three methods are presented and compared for Denton's 1-D nozzle. 

Calculation Details 

46 , <5x = 1 

1 , M = 0.80 

1.4 , R = 287. J/kg K 

0.85 , 0.80 , 0.75 


Figures 

The variations of static pressure, Mach number, and total pressure 
are plotted for each interpolation scheme using the same scales as for 
the theoretical solutions, Fig. 2. Fig. 4 shows the results for the 3- 


Number of Axial Grid Points = 
At inlet i = 
For air k = 

p exit/ p t inlet 


C -3- 


3.6 



point scheme, Fig. 5 for the 2-point scheme, and Fig. 6 for the M&M 
method. The results from the 3-point and M&M schemes are shown together 
with the theoretical solution on Fig. 7 for the pressure ratio of 0.80. 

Table 


The calculated values of maximum Mach number upstream of the shock 
and total pressure ratio are compared with the values from the theoreti- 
cal 1-D solutions in Table 2. 

Shock Losses 


The total pressure ratios across the shocks are well calculated by 
all three interpolation formulae as shown in Table 2b. This is in spite 
of the fact that the calculated values for the maximum Mach numbers 
upstream of the shocks are significantly different from the theoretical 
values. For example, at the lowest back pressure, the theoretical Mach 
number upstream of the shock is 1.578 while the 3-point interpolation 
formula gives 1.502, the 2-point formula gives 1.528, and the M&M for- 
mula gives 1.534. For this case the calculated values of total pressure 
ratio are all in the range 0.9027 to 0.9028, compared with the theoreti- 
cal value of 0.9032. In general the M&M formula gives the closest 
agreement with the upstream Mach number while the 3-point formula gives 
the worst results. Based on the maximum calculated upstream Mach number 
for these cases, the M&M formula would give shock losses from 16 to 42 
percent too small, while the 3-point formula would give values from 27 
to 62 percent too small. Interestingly the agreement for shock losses 
based on maximum upstream Mach number improves (for all three formulae) 
as the Mach number increases. However, these results show that peak 
calculated Mach number should not be used to predict shock losses and 
that the calculated total pressure loss across the shock is accurate to 
better than 0.1% and it should be used. 

Smoothing Upstream of Shock 

The results in Figs. 4, 5, and 6 show that the interpolation formu- 
lae all act to smooth properties upstream of the shocks. The smoothing 
is most noticeable in the static pressure and Mach number distributions, 
especially with the 3-point interpolation scheme. The 2-point scheme 
gives less smoothing while the M&M formula gives the sharpest and most 
accurate upstream distributions. 

Overshoots and Undershoots Downstream of Shock 


Both the 3-point and 2-point interpolation schemes give overshoots 
in static pressure and undershoots in Mach number downstream of the 
shocks. Only the M&M interpolation formula shows no noticeable over- 
shoots and undershoots and this is because it has a better formulation 
for subsonic flow; in fact, from Eq. 4, it can be seen that the M&M 
formula reduces to the correct perfect gas equation for Mach numbers 
less than 0.918. 
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Table 2. Results of Calculations for Denton's 1-D Nozzle. 
Table 2a: Maximum Mach number upstream of shock 

Interpolation Formula 


p 

exi t 

D 

Theoretical 

3-Point 

2-Point 

Mach Number 

t inlet 




Dependent 

0.85 

1.267 

1.173 

1.193 

1 .216 

0.80 

1.455 

1.375 

1.395 

1.408 

0.75 

1.578 

1.502 

1 .528 

1.534 


Table 2b: Total pressure ratio, P t exit /P t inlet 

Interpolation Formula 


P exit 

D 

Theoretical 

3-Point 

2-Point 

Mach Number 

t inlet 




Dependent 

0.85 

0.9847 

0.98487 

0.98492 

0.98494 

0.80 

0.9433 

0.94331 

0.94335 

0.94338 

0.75 

0.9032 

0.90271 

0.90277 

0.90281 
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Shock Location 


The M&M formula captures the shocks over about four grid points 
centered around the theoretical shock locations. This is seen for the 
pressure ratio of 0.8 in Fig. 7. In contrast, Fig. 7 shows the 3-point 
scheme smearing the shock over about ten grid points with the shock 
displaced slightly downstream due to inadequate resolution of the sub- 
sonic flow. Once again the 2-point scheme gives results intermediate 
between those of the M&M and 3-point schemes. 

3.5 Conclusions - Progress in Shock Capturing 

Significant improvements have been made in the finite-volume time- 
marching method to allow more accurate calculations of the distributions 
of flow properties through shocks. In Part 1, the Effective Density 
Method was introduced to reduce undershoots and overshoots in total 
pressure in the region of the shock. A stability analysis in Part 2 was 
then used to develop a Mach number dependent interpolation scheme for 
pressure which combines the advantages of the correct perfect gas equa- 
tion for subsonic flow with the stability of 2-point and 3-point inter- 
polation schemes for supersonic flow. The M&M interpolation formula, 
representing this new scheme, when used in the Effective Density Method, 
further removes the overshoot in static pressure in the subsonic flow 
downstream of a shock. 

References 

1. Denton, J. D. , "An Improved Time Marching Method for Turbomachinery 
Calculations,” ASME Paper 82-GT-239. 

2. Nicholson, S. , "Development of a Finite Volume Time Marching Meth- 
od," VPI&SU Turbomachinery Research Group Report No. JM/85-3, Febru- 
ary 1985. 

3. Nicholson, S., and J. Moore, "Extension of a Finite Volume Explicit 
Time Marching Method to Laminar and Turbulent Flow," VPI&SU Turbo- 
machinery Research Group Report No. JM/85-6, June 1985. 

4. Denton, J. D., "A Method of Calculating Fully Three Dimensional 
Inviscid Flow through Any Type of Turbomachine Blade Row," Lecture 
Notes for ASME Short Course on 3-D Flows in Turbomachinery Blade 
Rows, Phoenix, AZ, March 1983. 

5. Denton, J. D., "The Calculation of Fully Three Dimensional Flow 
through Any Type of Turbomachine Blade Row,” Lecture Notes for ASME 
Short Course on 3-D Flows in Turbomachinery Blade Rows, Houston, TX, 
March 1985, and AGARD Lecture Series No. 140 on 3-D Computation 
Techniques Applied to Internal Flows in Propulsion Systems, Rome, 
Italy, Cologne, West Germany, and Paris, France, June 1985. 


3.9 













MRCH NUMBER DEPENDENT R’S WITH R0+R1+R2=1 



0.00 0.20 0.40 0.60 0.80 1.00 1.20 l.MO 1.60 1.80 2.00 2.20 2.10 2.60 2.80 3.00 

MRCH 


Fig. 3 M & M Mach number dependent values (Eq . 4) 
for the coefficients in Eq. 2. 
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Fig. 4 Calculated 1-D solution for Denton's nozzle 
using 3-point interpolation, Eq. 2 with 
= 0 and a^ = 1. 
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Fig. 5 Calculated 1-D solution for Denton's nozzle 

using 2-point interpolation, Eq . 2 with 
a^ = a^ = 0 and a^ = 1. . 

Calculations for three exit static pressures at x = 46., 

P /p = 0.85, 0.80, and 0.75. 
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Fig. 6 Calculated 1-D solution for Denton's nozzle 
using M & M formula, Eq . 4. 

Calculations for three exit static pressures at x 

P . , = 0.85, 0.80, and 0.75. 
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Fig. 7 Comparison of calculated results with the theoretical 

1-D solution for P . ' /P . , =0.80. 
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